Esercitazione 4

Compensazione dinamica

Per la consegna dell’esercitazione si faccia riferimento a questo link. Mentre le librerie sono rimandate nell’icona Librerie.

1 Funzioni di utilità

# Generazione sinusoidale con diverse armoniche
signal <- function(t, pars, rad = FALSE) { 
  stopifnot(is.data.frame(pars))
  with(pars, {
    if (!rad) {
      phi <- phi/180*pi
      f <- 2*pi*f
    }
    map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*sin(t*f[i] + phi[i] ))))
  })
}

2 Simulazione uscita sistema di misura

# Parametri del misurando (input del sistema di misura)
pars <- tibble(
  w = c(1, 0.3, 0.3), 
  f = c(1/30, 1/10, 1/4), 
  phi = c(0, 0, 0)
)

# Generiamo il misurando u
Ts <- 0.01            # periodo di campionamento
fs <- 1/Ts            # freq. di campionamento   
u <- tibble(
  t = 0:10000 * Ts,
  u = signal(t, pars),
  un = u + rnorm(length(t), 0, pars$w[1]/10)
)

# Definiamo la Funzione di Trasferimento dello strumento
num <- c(5, 1)   # 5s + 1
den <- c(1, 2, 1) # s^2 + 2s + 1
H <- tf(num, den)
print(H)

y1:
    5 s^1 + 1 
  - - - - - - - - - -
     s^2 + 2 s + 1 


Transfer Function: Continuous time model 

# Simulazione dell'uscita nel tempo
tmp <- u %>% mutate(
    out = lsim(H, un, t)$y[1,],
    ideal = lsim(H, u, t)$y[1,]
  ) %>% mutate(n = row_number(), .before = t)

tmp %>% plot_ly(x = ~t) %>%
  add_lines(y = ~un, name = "Input") %>%
  add_lines(y = ~ideal, name = "Output ideale") %>%
  add_lines(y = ~out, name = "Output") %>%
  layout(xaxis=list(title="Tempo (s)"), yaxis=list(title=""), title = "Simulazione dell'uscita")

Notare che lsim calcola la risposta nel tempo di un sistema lineare.

2.1 Spettro dei segnali misurati

Ta <- max(u$t)
un.fft <- compute_fft(tmp, Ta, un)
ideal.fft <- compute_fft(tmp, Ta, ideal)
out.fft <- compute_fft(tmp, Ta, out)

bind_rows(
  un.fft %>% mutate(caso = "Input"),
  ideal.fft %>% mutate(caso = "Output ideale"),
  out.fft %>% mutate(caso = "Output")
) %>% dplyr::filter(f < 1) %>% 
  ggplot(aes(x = f, y = mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#F8766D") +
  geom_point(colour = "#F8766D") +
  facet_wrap(~ caso, ncol = 1, scales = "free_x") +
  labs(x = "Frequenza (Hz)", y = "Intensità")

NB: non si vede, ma il rumore bianco è presente. Provare ad ingrandire la sd del rumore addizionato a un.

3 Stima del misurando a partire dal segnale di uscita

# Segnale in uscita con rumore
yI <- lsim(H, u$u, u$t)$y[1,]
yn <- yI + rnorm(length(u$t), 0, pars$w[1]/10)


# Ricaviamo l'inversa della Funzione di Trasferimento:
Hinv <- tf(den, num)

[1] "TFCHK: Transfer function may not be proper and may lead to errors. Num > Den"

print(Hinv)

y1:
     s^2 + 2 s + 1 
  - - - - - - - - - -
    5 s^1 + 1 

Transfer Function: Continuous time model 

# Simulazione dell'ingresso
#input <- lsim(Hinv, yn, y$t, x0 = rep(0, length(pole(Hinv))))

Error in tf2ss(num, den) :
TF2SS: The order of the Numerator should be equal or lesser than the Denominator.

Implementazione della compensazione dinamica.

Il problema risiede nel fatto che la funzione lsim è progettata per calcolare la risposta nel tempo di un sistema LTI proprio (\(deg(\text{NUM}) \leq deg(\text{DEN})\)). Nel caso dell’inversa \(H^{-1}(j\omega)\), è come se si stimasse la risposta di un sistema improprio, per cui lsim fallisce.

La soluzione impiegata è quella di calcolare numericamente il misurando nel dominio della frequenza, dopodiché mediante l’anti-trasformata si trova il rispettivo segnale nel tempo. Pertanto, i passaggi da effetuare saranno:

1. trasformare l’uscita nel dominio della frequenza;
2. calcolare la risposta in frequenza della FdT;
3. calcolare il misurando invertendo la caratteristica dinamica: \(U(j\omega) = \frac{Y(j\omega)}{H(j\omega)}\);
4. anti-trasformare per ottenere \(u(t)\).

# Filtraggio dell'uscita per ridurre il rumore
y <- u %>% select(t) %>% mutate(n = row_number(t), .before = t) %>%
  mutate(
    y_id = yI,
    y_n = yn
  )

tmp <- compute_fft(y, Ta, yn)

tmp %>% dplyr::filter(f < 2.5) %>%
  ggplot(aes(x = f, y = mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#56B4E9") +
  geom_point(colour = "#56B4E9") +
  labs(x = "Frequenza (Hz)", y = "Intensità", title = "Spettro dell'uscita con rumore")

# Filtro Butterworth con taglio a 1Hz
fc <- 1
flt <- butter(2, fc/(fs/2))

#ggbodeplot_digital(flt, fs, fmin = 1e-1, fmax = 1e1) + geom_vline(xintercept = fc, color="red", linetype = 2)

y <- y %>% mutate(y_flt = signal::filtfilt(flt, yn))

# `filtfilt` effettua di default un padding ai bordi del segnale,
# introducendo una distorsione nella coda destra.
# Perciò viene tagliata tale coda
y <- y %>% dplyr::filter(t <= 99)

# Nuovo tempo di campionamento
Ta <- max(y$t)

Si osservi come la banda utile dell’uscita si sovrapponga alla banda di stabilità della FdT inversa. Questo darà la garanzia che il rumore ad alta frequenza non venga amplificato.

N <- length(y$t)

# Asse frequenze
freq <- tibble(
  k = 0:(N-1),
  f = ifelse(k <= N/2, k, k-N) * fs/N,   # FFT simmetrica all'asse delle freq.
  w = 2*pi*f
)

# Trasformata di Fourier (FFT)
y <- y %>% mutate(Y =fft(y_flt))

# Calcolo della funz. di risposta in frequenza
Hiw <- freqresp(H, freq$w)

# Inversione della caratteristica dinamica e anti-trasformata
u_comp <- tibble(
  t = y$t,
  U = y$Y / Hiw,
  u = Re(fft(U, inverse = TRUE)) / N  # Re() serve a trascurare dei residui di calcolo nella parte Im. 
                                      # La fft(inverse = TRUE) ritorna l'inversa NON-normalizzata, 
                                      # perciò si divide per la lunghezza del campione
)

Il risultato ottenuto è molto buono nella parte centrale. Tuttavia, all’inizio si osserva una discrepanza tra i dati compensati e quelli desiderati e delle forti oscillazioni su entrambi i bordi. La ragione risiede nel fatto che la FFT assume che il segnale sia periodico e campionato per un numero di cicli interi (Finestre), evitando discontinuità e non derivabilità ai bordi. Tali condizioni non sono soddisfatte in questo caso.

Pertanto, si possono implementare due possibili soluzioni:

• [proposta da ChatGPT] effetuare un’operazione di finestratura del segnale;
• effettuare un padding a zero sulle code che impone la condizione \(u(0) = u(T_a) = 0\).

Prima si prova ad effettuare il windowing.

# Windowing
y <- y %>% mutate(
  win = signal::hanning(N),
  y_win = y_flt * win
)

# Trasformata di Fourier (FFT)
y <- y %>% mutate(Y_win = fft(y_win))

# Calcolo della funz. di risposta in frequenza
Hiw <- freqresp(H, freq$w)

# Inversione della caratteristica dinamica e anti-trasformata
u_comp <- u_comp %>% mutate(
  U_win = y$Y_win / Hiw,
  u_win = Re(fft(U_win, inverse = TRUE)) / N
)

Si può concludere che l’operazione di finestratura non è la più ottimale: la ricostruzione del misurando risulta buona solo nel centro della finestra, risultato che era stato già ottenuto precedentemente, anzi anche migliore. Il fatto è che, nel caso di questa applicazione, si preferisce ricostruire l’intero segnale di ingresso e non solo una parte.

3.1 Padding simmetrico a zero

# Ampiezza padding
M <- 2*N

pad_left  <- floor((M - N)/2)
pad_right <- ceiling((M - N)/2)
t_start <- y$t[1] - pad_left * Ts

y_ext <- tibble(
  t = seq(from = t_start, by = Ts, length.out = M),
  y_pad = c(
    rep(0, pad_left),
    y$y_flt,
    rep(0, pad_right)
  )
)

# Verifica che il segnale sia centrato: deve restituire TRUE
#all.equal(y_ext$y_pad[(pad_left+1):(pad_left+N)],y$y_flt)

# Nuovo asse delle frequenze
freq <- tibble(
  k = 0:(M-1),
  f = ifelse(k <= M/2, k, k-M) * fs/M,   # FFT simmetrica nell'asse delle freq.
  w = 2*pi*f
)

Hiw <- freqresp(H, freq$w)

y_ext <- y_ext %>% mutate(
  Y_pad = fft(y_pad),
  U_pad = Y_pad/Hiw,
  u_pad = Re(fft(U_pad, inverse = TRUE)) / M
)

# Scarto delle code di padding
u_comp <- u_comp %>% mutate(
  u_pad = y_ext$u_pad[(pad_left+1):(pad_left+N)]
)

Questa volta il segnale compensato coincide molto bene con quello ideale. Ciononostante, si osserva ancora un alto contenuto in frequenza nelle code. Quest’ultimo viene eliminato mediante un’operazione di filtraggio.

# Verifico lo spettro del segnale
compute_fft(u_comp, Ta, u_pad) %>% dplyr::filter(f < 2) %>%
  ggplot(aes(x = f, y = mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#56B4E9") +
  geom_point(colour = "#56B4E9") +
  labs(x = "Frequenza (Hz)", y = "Intensità", title = "Spettro dell'ingresso compensato")

fc <- 2
flt <- butter(2, 2*fc/fs)

u_comp <- u_comp %>% mutate(
  u_flt = signal::filtfilt(flt, u_pad)
)

Adesso si osserva che l’ingresso viene riprodotto fedelmente. Di seguito un confronto tra segnale di uscita e ingresso compensato: